In this lab you will learn



library(tidyverse)
library(wordspace) ## for function dist.matrix()
library(ggforce) ## for functions geom_circle() and coord_fixed()
theme_set(theme_bw())
theme_update(
  aspect.ratio = 1,
  panel.grid = element_blank()
)



Metapopulation theory aside, some facts remain: the average species has many subpopulations, many species often use discontinuous habitat, and dispersal seems to influence population viability. Shafer 2001, Biological Conservation

In this lab, we will examine how the population dynamics of a species can be affected when it lives in disjointed patches (or “islands”) of suitable habitat rather than in a single contiguous region. Specifically, we are interested in the metapopulation scenario, where each patch or island may or may not contain a reproductive population at any given time, individuals can occasionally disperse between islands (immigration), and local extinction and recolonization events are possible for all islands.

Keep in mind that the term island in the context of island biogeography theory and metapopulation theory is broadly interpreted as a contiguous patch of suitable habitat surrounded by unsuitable (or less suitable) habitat, such that there is occasional exchange of individuals between these populations. Therefore, “islands” here can mean mountain tops, coral reefs, streams, ponds, parks, and even bodies–as in hosts for a pathogen!


The theoretician’s perspective of metapopulations, by Rampal S Etienne, doctoral dissertation


Thus broadly defined, many species in nature comprise metapopulations, especially as result of anthropogenic habitat fragmentation. Understanding metapopulation dynamics is therefore of extreme importance to conservation and management.



Figure: Examples of metapopulations, defined as populations distributed over disjointed or fragmented habitat. The western fence lizard (Sceloporus occidentalis) has two subpopulations in Washington State (green/yellow patches). The American pika (Ochotona princeps) occurs in multiple subpopulations distributed across the American West. The holm oak (Quercus ilex) is found along the Mediterranean coast, plus isolated populations across Western Europe and Turkey. Photo credits by fksr (Flickr), Nature Mapping Foundation, Frédéric Dulude-de Broin, Chermundy, Gardenia.net, Giovanni Caudullo.



Metapopulation model

Suppose we have an archipelago with 20 islands spanning approximately 1,000 km by 1,000 km, roughly similar to the area spanned by the Hawaiian Windward Islands. For simplicity we’ll assume all islands are round, and for now let’s give them all equal diameter of 30km, as wide as Long Island and about the size of Maui.

Figure: Left: Hawaiian Windward Islands, by Jacques Descloitres. Right: The islands in our model.

Now, let’s think about how our model should work. Each island contains suitable habitat for our focal species – say, a finch. Each year, the local population on an occupied island can go extinct through chance events – a food shortage, or a bout of disease or catastrophic weather. Also each year, an island which is currently unoccupied by the focal species may be recolonized via wayward immigrants from occupied islands. This is a verbal description of Levins’s metapopulation model,

\[ \frac{dp}{dt} = cp(1-p) - ep \]

where \(p\) is the occupancy (proportion of islands currently occupied), and \(c\) and \(e\) are the rates of immigration and extinction, respectively. Our model will take into account the factors that affect these rates, and will also take into account the fact that nature is probabilistic rather than deterministic as in the equation above.


Q1. Based on the definition of occupancy given above, if an archipelago has 20 islands and the current occupancy is 0.4, in how many islands is the species currently present?


What factors should affect the yearly probability of local extinction on an island? One of the most important factors is population size: the larger the population, the lower the chance of extinction. Population size in turn depends on habitat size and habitat quality. We will assume all islands are equally suitable for the species. We will thus write our extinction term as an extinction constant reflecting habitat quality (the same for all islands since we’re assuming they are all equally suitable), divided by the island’s area.

What factors should affect the yearly probability of a currently vacant island being recolonized by immigrants from other islands? There are four main factors to consider:

So our immigration term is the product of an intrinsic immigration rate, the area of the source island, the area of the target island, and the square inverse of the distance between them, summed over all possible source islands (those currently occupied).

The theoretical literature on metapopulations provides predictions for the probability of persistence of a metapopulation given the factors above (Hanski & Ovaskainen 2000, Grilli et al 2015). However, here we will follow an empirical approach, where we numerically simulate a metapopulation model and check how persistence is affected by our parameters of interest (habitat size, connectivity, etc).

Here is the R code for our model:

Metapop_Model = 
  function(
    number_of_islands,
    final_year,
    extinction_constant,
    immigration_constant,
    areas,
    coordinates,
    initial_state = rep(1, number_of_islands),
    seed = 0
  ){
    # Standardize the random sequence of events
    set.seed(seed)

    # Derived parameters
    distances2 = as.matrix(dist(coordinates) ^ 2)
    extinction_rates = extinction_constant / areas
    immigration_flux = areas / distances2; diag(immigration_flux) = 0
    
    # Initial conditions
    extinction_year = NA
    state = initial_state
    occupancy = mean(state)
    immigration_rates = immigration_constant * areas * colSums((state == 1) * immigration_flux)
    state_tibble = 
      enframe(state) |> 
      pivot_wider(
        names_from = name, 
        values_from = value
      )
    colnames(state_tibble) = paste0('island_', seq(number_of_islands))
    
    # Metapopulation dynamics
    for(year in 1:final_year){
      
      extinction = rbernoulli(state, p = extinction_rates * (state == 1))
      state[extinction] = 0
      
      immigration = rbernoulli(state, p = immigration_rates * (state == 0))
      state[immigration] = 1
      
      immigration_rates = immigration_constant * areas * colSums((state == 1) * immigration_flux)
      
      occupancy = c(occupancy, mean(state))
      
      state_tibble = rbind(state_tibble, state)
      
      if(all(state == 0)){
        extinction_year = year
        break
      } 
    }
    
    time_at_end_of_simulation = ifelse(is.na(extinction_year), final_year, extinction_year)
    
    # Measurements
    state_tibble = 
      bind_cols(
        tibble(year = 0:time_at_end_of_simulation),
        state_tibble
      )
    
    time_averaged_occupancy = mean(occupancy[-(1:100)])
    
    return(
      list(
        parameters = 
          list(
            number_of_islands = number_of_islands,
            final_year = final_year,
            e = extinction_constant,
            c = immigration_constant,
            coordinates = coords,
            areas = areas
          ),
        metapopulation = state_tibble,
        occupancy = occupancy,
        global_extinction = all(state == 0),
        time_at_end_of_simulation = time_at_end_of_simulation
      )
    )
  }


Time-averaged occupancy

Let’s introduce R code that takes the results of the model and plots the occupancy over time:

Plot_Occupancy = function(model){
  p = model$occupancy
  time_at_end_of_simulation = model$time_at_end_of_simulation
  time_averaged_occupancy = round(mean(p[-seq(time_at_end_of_simulation / 2)]), 2)
  
  plot = 
    tibble(
      year = seq_along(p),
      occupancy = p
    ) |>
    ggplot(aes(year, occupancy)) +
    geom_line() +
    coord_cartesian(ylim = c(0, 1)) +
    ggtitle(
      bquote(
        paste(
        'number of islands = ', .(model$parameters$number_of_islands),
        '    e = ', .(model$parameters$e),
        '    c = ', .(model$parameters$c),
        '    ', bar(p), ' = ', .(time_averaged_occupancy)
        )
      )
    ) +
    geom_hline(aes(yintercept = time_averaged_occupancy), color = 'red') +
    geom_vline(aes(xintercept = time_at_end_of_simulation / 2), color = 'blue', linetype = 'dotted')
    
  plot |> 
    show()
}


Now we must parametrize our model. Let’s set the extinction constant to 100 and the immigration constant to 0.0015. Given our other model parameters (island sizes and distances), this works out to an average annual probability of extinction of 14% per occupied island, and average annual probability of recolonization of 10% per unoccupied island (the latter will vary depending on the current set of occupied islands).

# standardize the random number generator
set.seed(30)

# island parameters
number_of_islands = 20
areas = rep(pi * 15 ^ 2, number_of_islands)
coords = 
  tibble(
    x = runif(number_of_islands, min = 0, max = 1e3),
    y = runif(number_of_islands, min = 0, max = 1e3)
  )

# immigration parameter
immigration_constant = 0.001

# extinction parameter
extinction_constant = 100


Q2. What ecological factors can affect the extinction rate of a local population within a patch of habitat of a given size?


Q3. What ecological factors can affect the immigration rate into an island of a given size and distance to other islands?

We are now ready to run the model. Let’s imagine all islands are occupied in year 0, and run the model until year 100. Here’s what the occupancy over time looks like:

model = 
  Metapop_Model(
    number_of_islands = number_of_islands,
    final_year = 100,
    extinction_constant = extinction_constant,
    immigration_constant = immigration_constant,
    areas = areas,
    coordinates = coords,
    seed = 0
  )

Plot_Occupancy(model)

This population persisted through the end of the simulation, and remains alive after 100 years!

The red line shows the average occupancy taken across the latter half of the years for which we ran the simulation (\(\bar{p} = 0.37\)). The halfway mark is shown by the dotted blue line.


Q4. Why not take the average occupancy all the way from the beginning of the simulation at year 0?

For Q5 and Q6, use the code provided below, replacing the ?? with the appropriate parameter values.

Plot_Occupancy(
  Metapop_Model(
    number_of_islands = 20,
    final_year = 100,
    extinction_constant = ??,
    immigration_constant = ??,
    areas = areas,
    coordinates = coords,
    seed = 0
  )
)


Q5. What happens to the time-averaged occupancy when we double the immigration constant from 0.001 to 0.002? Rerun the code above after changing the appropriate parameter. Show a plot similar to the one above for this scenario.


Q6. What happens when we keep the immigration constant at 0.0015 but reduce the extinction constant from 100 to 50? Show a plot similar to the one above for this scenario.


Q7. What happens when we destroy the habitat in half of the islands? Does the population persist for 100 years? Find out by running the code below, which simulates the model again, this time using the original values of both the extinction and immigration constants but randomly selecting 10 of the original 20 islands as remaining viable habitat for our species. Show a plot similar to the one above for this scenario.

Plot_Occupancy(
  Metapop_Model(
    number_of_islands = 10,
    final_year = 100,
    extinction_constant = 100,
    immigration_constant = 0.001,
    areas = areas[1:10],
    coordinates = coords[1:10, ],
    seed = 0
  )
)


You may have noticed from Q5 and Q6 that the per-island rates of extinction and immigration greatly affect the proportion of habitat we expect to be occupied by our species over time. From a modeling perspective this makes sense, as those influences were directly built into the model. When we doubled immigration, we got twice the average occupancy, which may seem reasonable.

However, notice the drastic difference it made to remove habitat in Q7! By destroying the habitat in half of the islands, we lost the entire metapopulation almost immediately! This happened even though less than half of the original habitat was usually occupied in the first place.


Q8. One of the advantages of ecological modeling is that we can learn new things by examining unanticipated results with the benefit of hindsight. Considering the way we built the model, what are some of the reasons why habitat loss had such an outsize influence on the outcome for our species?


Another important factor affecting the fate of the metapopulation is stochasticity. As you may have noticed, this model contains a random element to it, in the sense that immigration and extinction are probabilistic events. This means that if we rerun the model with the same parameters and initial conditions, we may get a different result.


Q9. Let’s check this statement by rerunning the model, this time with a different seed for the random number generator. Set the seed argument in the call to function Metapop_Model to 10, while setting all other parameters back to the original scenario. Do you see the same occupancy curve as in the plot above?


You may now be wondering what would happen if we repeated the experiment many times. How variable will the average occupancy be? Is the mean occupancy we got the first time representative of what is most likely to happen? Might we even get global extinctions in some of the runs?

Let’s check. Copy/paste the code chunk below, which defines a function which runs the model many times and plots the results.

Multiple_Sims = 
  \(
    .numruns,
    .final_year,
    .number_of_islands,
    .extinction_constant,
    .immigration_constant,
    .areas,
    .coordinates
  ){
    results = NULL
    for(seed in 1:.numruns){
      model = 
        Metapop_Model(
          number_of_islands = .number_of_islands,
          final_year = .final_year,
          extinction_constant = .extinction_constant,
          immigration_constant = .immigration_constant,
          areas = .areas,
          coordinates = .coordinates,
          seed = seed
        )
      
      results = 
        results |>
        bind_rows(
          tibble(
            seed = seed,
            time_at_end_of_simulation = model$time_at_end_of_simulation,
            global_extinction = model$global_extinction,
            time = seq(0:.final_year),
            occupancy = c(model$occupancy, rep(0, 1 + .final_year - length(model$occupancy)))
          )
        )
    }
    
    extinctions = 
      results |> 
      select(seed, global_extinction) |> 
      unique() |> 
      summarize(
        extinctions = sum(global_extinction), 
        .groups = 'drop'
      ) 
    
    plot = 
      results |> 
      group_by(time) |> 
      summarize(
        p = mean(occupancy), 
        sd = sd(occupancy),
        .groups = 'drop'
      ) |> 
      ggplot(aes(time, p)) + 
      geom_ribbon(aes(ymin = p - sd, ymax = p + sd), fill = rgb(87.8/100, 69/100, 1, .5)) +
      geom_line() +
      geom_text(
        aes(
          x = .75 * .final_year, 
          y = .9, 
          label = paste(extinctions, ' global extinctions')
        ), 
        data = extinctions
      ) +
      labs(x = 'Year', y = 'Occupancy') +
      ggtitle(
        bquote(
            paste(
            'no. islands = ', .(model$parameters$number_of_islands),
            '    e = ', .(model$parameters$e),
            '    c = ', .(model$parameters$c)
            )
          )
      )
    
    return(list(results = results, plot = plot))
  }


The plot below was generated by running the model 100 times, each with a different seed (i.e. with a different set of random numbers generated during the simulation). For each year, it shows the average occupancy across all simulations for that year (black line). The lilac band around the line gives a measure of the variation around the mean.

runs = 
  Multiple_Sims(
    .numruns = 100,
    .final_year = 500,
    .number_of_islands = 20,
    .extinction_constant = 100,
    .immigration_constant = .001,
    .areas = areas,
    .coordinates = coords
  )

runs |>
  pluck('plot') |>
  show()

We notice that mean occupancy immediately drops from 100% to just about 25% in the first 50 years of the simulations, then continues to drop at a slower rate. Notice the high variability around the mean, which indicates that in some simulations the occupancy at a given point in time is much higher than in others.By the 500th year, the entire metapopulation has gone extinct in 73 out of the 100 runs.


Q10. What do you think would happen in the 27 simulations where the population was still alive by Year 500 if we run kept running them beyond 500 years? (No need to run any code to answer this question.)


Let’s now repeat the experiment above, this time with an immigration constant of 0.0005, which is one half of the original value we used.


Q11. Perform the task above by copy/pasting the code chunk above the plot above, and changing the value of the immigration constant to 0.0005. Also, change the .final_year argument to 200. Compare your plot to the plot above and explain the differences in as much detail as you can.

The histogram below plots the distribution of the number of years that it took for a simulation to reach global extinction in the scenario your just ran in Q11:

There is a big spread in global extinction times: in some runs the species lasts for 100 years, in others it gets wiped out in less than 25 years. This happened despite the fact that every parameter in the model had the same value across all runs. Why this huge variation? Because when the total population of the species is low enough to make it highly vulnerable to extinction, stochastic forces play an outsize role, and the system becomes more unpredictable.

The lesson to draw here is profound: because drift creates a non-zero probability that the entire species will go extinct from one year to the next, then we can say with certainty that the species will go extinct — even if it takes a long time. Notice that this applies to all species, including our own: since extinction of the human species is not impossible, then it is certain to happen at some point in the future!

Figure: Sherlock Holmes’s famous quote on stochastic processes.


Now that we convinced ourselves that every species ultimately goes extinct, the question is how long until it happens, or how likely it will happen within a specified time window. This is why in conservation biology the better question is not “Is species X vulnerable to extinction”, but “What is the probability that species X will go extinct within Y years under current circumstances?” If we have a good understanding of how the circumstances affect this extinction risk, we can try to change those circumstances to bring the extinction risk down to acceptable levels.


Extinction risk

Let’s examine the model to see how the probability of extinction depends on factors that policy makers may be able to do something about. Here we will focus on habitat size and island connectivity (again, island here is meant sensu latu.)

First let’s take a closer look at the occupancy of each island over time. Let’s run the model with the original parameters for 1,000 years, and take note of what percentage of time each island was occupied.

model = 
  Metapop_Model(
    number_of_islands = 20,
    final_year = 1000,
    extinction_constant = 100,
    immigration_constant = 0.0015,
    areas = areas,
    coordinates = coords,
    seed = 0
  )

island_occupancy_tibble =
  bind_cols(
    coords,
    island_occupancy = as.numeric(round(100 * colMeans(model$metapopulation[, -1])))
  ) |>
  rowid_to_column(var = 'island_index')

plot_occupancy_by_island = 
  island_occupancy_tibble |>
  ggplot() +
  geom_point(aes(x, y, color = island_occupancy), alpha = .7, size = 8) +
  geom_text(aes(x, y, label = island_occupancy)) +
  scale_color_gradient2(midpoint = 50, low = "green", mid = 'yellow', high = "red") +
  theme(
    legend.position = 'none',
    aspect.ratio = 1,
    panel.background = element_rect(fill = 'lightblue')
  )


The map below shows our islands again, this time colored and labeled by the percentage of years where they were occupied by our focal species.

plot_occupancy_by_island

Q12. Are all islands populated equally often? What pattern do you see? How do you explain the pattern?


Q13. Suppose the islands in our archipelago are suitable for growing a certain commercial crop, and policymakers decided to allow clearing the native habitat of three of the islands for agriculture. Aiming to minimize the impact on our focal species, the local environmental agency hires you as a consultant. a) Which three islands would you designate for clearing? b) Conversely, which three islands would you strongly recommend against clearing? Check your answer by rerunning the model above without the islands on your best- and worst-candidates lists and see what happens, which you can do using the code provided below. To find out which islands you’d like to exclude, look up the island_index column of the island_occupancy_tibble object, and substitute those indices for the ?? in the code below (for example, if the best islands for removal – i.e. those whose loss will have the least effect on our species – are islands 4, 13, and 12, then set index_1 = 4, index_2 = 13, and index_3 = 12, then run the code). Show your plots for (a) and (b).

index_1 = ??
index_2 = ??
index_3 = ??
  
Plot_Occupancy(
  model = 
    Metapop_Model(
      number_of_islands = 17,
      final_year = 1000,
      extinction_constant = 100,
      immigration_constant = 0.0015,
      areas = areas[-c(index_1, index_2, index_3)],
      coordinates = coords[-c(index_1, index_2, index_3), ],
      seed = 0
    )
)


Figure: Oil palm plantation in Papua New Guinea. Photo by Jurgen Freund, WWF-Canon


Note: Other things being equal, this experiment gives us a sense of the importance of habitat connectivity for metapopulations. Of course, things tend to be much more complicated in real life, where other things are seldom equal. For example, we are assuming all islands are equally suitable habitat. Even if that is true for our focal species, it may not be true for other species.


The SLOSS debate

Conservation biologists debated in the 1970s and 1980s whether a single large reserve is better than several small ones with equivalent total area. Jared Diamond suggested that one large protected area is more likely to protect species from extinction, whereas Daniel Simberloff and others countered that in some circumstances, several small reserves could be more effective – for example, small islands spanning a wider spatial range could contain a higher variety of microhabitats, affording species a larger number of local niches. This became known as the SLOSS debate (for Single Large Or Several Small).

Figure: Alternative reserve designs for a total given area. Credit: Brittany Mendrikis.


Q14. In Q12, we considered habitat islands separated from each other by water. Now imagine what would happen if we took all of the islands and amalgamated them into one single island with total area equal to the sum of the areas of the smaller islands. The total habitat area would not change, but the average connectivity would. a) Would connectivity in this single island increase or decrease, relative to the average connectivity among the islands in Q12? b) Considering what you saw in Q12 regarding the impact of connectivity on species persistence, would you expect higher species persistence in this single island or in the archipelago in Q12? c) What side of the SLOSS debate is supported by these ideas?


Now let’s relax the assumption that all islands have the same area. Consider a scenario where our 20 islands differ in size, while averaging the same area between them as our simpler identical-areas scenario.


Q15: Which scenario (all islands have same size versus islands have different sizes) do you predict would give our species a better chance at long-term survival? Explain your reasoning.

The code below will test your prediction. The code runs the identical-areas scenario for 500 years with an immigration constant of 0.0005, and records the time to global extinction in each of 100 runs. Then it repeats the experiment under a scenario where island areas are distributed like a bell curve (also called the normal distribution) with a standard deviation of 40% of the mean. Finally, it plots the distribution of times to extinction in each scenario as a density plot, which is basically a histogram on a continuous axis.


set.seed(10)

result_identical_areas = NULL

for(run in 1:100){
  model = 
    Metapop_Model(
      number_of_islands = 20,
      final_year = 500,
      extinction_constant = 100,
      immigration_constant = 0.0005,
      areas = areas,
      coordinates = coords,
      seed = run
    )
  
  result_identical_areas = 
    c(result_identical_areas, model$time_at_end_of_simulation)

}


result_different_areas = NULL

for(run in 1:100){
  model = 
    Metapop_Model(
      number_of_islands = 20,
      final_year = 500,
      extinction_constant = 100,
      immigration_constant = 0.0005,
      areas = pmax(pi * (15/20) ^ 2, rnorm(n = 20, mean = pi * 15 ^ 2, sd = .5 * pi * 15 ^ 2)),
      coordinates = coords,
      seed = run
    )
  
  result_different_areas = c(result_different_areas, model$time_at_end_of_simulation)

}

plot_results = 
  tibble(
    identical_areas = result_identical_areas,
    different_areas = result_different_areas
  ) |>
  rowid_to_column(var = 'run') |>
  pivot_longer(-run, names_to = 'scenario', values_to = 'time_to_extinction') |>
  ggplot(aes(time_to_extinction, group = scenario, fill = scenario, color = scenario)) +
  geom_density(alpha = .3)

plot_results


Q16: Which scenario is more variable (i.e. where do you see a wider spread of times to extinction)? Under which scenario is the species more likely to still be alive after 500 years? Could you have foreseen these results without the model?


Reserve design

The exercises above gave you a sense of what managers and conservationists must consider when planning wildlife reserves. Of course, our model is very simple. In reality, modelers must attempt to assess and incorporate the complexities germane to the species one wishes to protect. Even then, one must always contend with issues that have little to do with metapopulation dynamics, and yet are equally important for conservation success.


Q17. In addition to habitat quality and connectivity, name one factor that may affect the prospects of long-term survival of a metapopulation.


One factor to consider is that the greatest risks to persistence for some species may not relate to reserve design. For example, population viability analyses of the mountain gorilla (Gorilla gorilla beringei) estimate that the species is unlikely to go extinct solely due to genetic or demographic reasons (factors affected by reserve size, number, and connectivity) for several hundred years. Rather, the main threat to gorillas is poaching, a problem whose resolution lies in patrolling existing reserves and conservation education.

Figure: Mountain gorilla and Sumatran rhinoceros, two species whose future may rely on poaching control.


Despite being extrinsic to metapopulation dynamics sensu strictu, the threat from poaching leads to yet another concern in reserve design. While larger and more numerous reserves protect against stochastic extinction and genetic inbreeding, they are also more difficult to patrol against poaching. For example, “the endangered Sumatran rhinoceros (Dicerorhinus sumatrensis) persists at about 50 sites in Asia but such widespread distribution has been judged too difficult to protect against poaching” (Shafer 2001).

On the other hand, distance between reserves may protect subpopulations from catastrophic events, such as hurricanes and epidemics. For example, in 1989 Hurricane Hugo destroyed 63% of the local population of the red-cockaded woodpecker (Dryobates borealis) in South Carolina’s Francis Marion National Forest.

Figure: The red-cockaded woodpecker (Dryobates borealis), once distributed across the entire southeastern US and beyond, is currently found in fragmented patches.


Q18. The red-cockaded woodpecker is currently distributed across the American Southeast in fragmented habitat, and is classified as Near Threatened by the IUCN. What could have happened if its entire remaining population had been concentrated near South Carolina prior to 1989?

While distance between local populations may protect against localized catastrophes, we have seen in this lab that habitat connectivity is crucial to support a metapopulation.

One mitigation effort against separation between reserves is the construction of wildlife corridors. Corridors may allow animals to escape fire and poaching, in addition to the benefits from high connectivity reviewed in the exercises above. Unfortunately, it is not easy to measure the effectiveness of corridors in connecting populations, and the available evidence is mixed. This is not per se an argument against corridors (and there are many arguments for them, see eg Noss 1987), but it further highlights that reserve design must contend with limited data.

Furthermore, habitat connectivity also has its downsides.

From Rampal S Etienne, doctoral dissertation


Q19. Based on the cartoon above, explain one disadvantage of reserve connectivity.

—